//lanczos.cxx

#include <iostream>
#include <cmath>
#include "lanczos.h"
#include "myvectormath.h"
#include "eivtris.h"
#include "global.h"
using namespace std;

int lanczos(Int_t N, vector<Double_t> *s, void (*matvec)(Int_t n, vector<Double_t> *u, vector<Double_t> *v)){
  
  vector<Double_t> q(N); 
  vector<Double_t> v1(N); // Verschiebungsvektor
  vector<Double_t> v2(N); // multi purpose temporay vector
  vector<Double_t> alpha;
  vector<Double_t> beta;
  vector<Double_t> temp1(N);
  vector<Double_t> temp2(N);
  
  Double_t normnew;
  
  normnew=skalarp(s,s);
  if (normnew==0){
    cout << "Norm des Startvektors ist 0." << endl;
    return 0;
  }
  VectorMultiS(s, 1/normnew, &v1);
  matvec(N,&v1,&q);
  
  for (int i=0; i<N;i++){
    alpha.push_back(skalarp(&v1,&q));
    VectorMultiS(&v1,alpha.at(i),&temp1);
    vectorsub(&q,&temp1,&q);
    beta.push_back(skalarp(&q,&q));
    if (beta.at(i)==0){
      cout << "beta(" << i << ")=0" << endl;
      break;
    }
    VectorMultiS(&q,beta.at(i),&v2);
    matvec(N,&v2,&temp1);
    VectorMultiS(&v1,beta.at(i),&temp2);
    vectorsub(&temp1,&temp2,&q);
  }
  
//   int* array = new int[length];
  Double_t * alphaArray = new Double_t[N];
  Double_t * betaArray = new Double_t[N-1];
  Double_t * eival = new Double_t[N];
  
  eivtris(N, alphaArray, betaArray, eival);
  cout << "calculated eigenvalues: \n";
  for (int i=0; i<N;i++){
    cout << eival[i] << " ";
  }
  cout << endl;
  
  delete alphaArray;
  delete betaArray;
  delete eival;
  
  return 0;
}